% This file take the steady-state the file "toDynare.mat"
% It writes the dynare file "code_dynare_hat.mod" which is the dynamic model
% to use perturbation techniques
% it double-checks that the steady state is OK, thanks to the "resid" function




clear
warning('off','all');
isOctave = exist('OCTAVE_VERSION', 'builtin') ~= 0;

isOctave && warning('off', 'Octave:array-as-logical');

load toDynare

trunc  = Ramsey.truncatedModel.truncatedAllocation;
lagMul = Ramsey.lagrangeMult;

Ntot  = Ramsey.truncatedModel.Ntot;    %total number of bins
Pi_h  = trunc.Pi_h';
S_h   = trunc.S_h;
K     = solution.K(1);
Ltot  = solution.L(1);
alpha = eco.alpha(1);
Y     = K^alpha*Ltot^(1-alpha);
FL    = (1-alpha)*K^alpha*Ltot^(-alpha);
xsis  = Ramsey.truncatedModel.xsis;

%xsis.xsiu1 = xsis.xsiuE;
xsis.xsiu2 = xsis.xsiuE;



ys    = eco.ys;


rho_u = 0.1;
G0 = .01*(1-rho_u); 


l     = Ramsey.l;
taut =  (( 1.0/eco.phi + 1.0)*eco.tau)/( 1.0/eco.phi +1.0 - eco.tau);
x = trunc.c_h - (1/(eco.chi*(1/eco.phi + 1)))*(trunc.l_h).^(1+1/eco.phi);


P     = ones(Ntot,1);
P(trunc.ind_cc_h) = 0;


file = 'code_dynare.mod';
fid  = fopen(file, 'w');
formatSpec = '%25.20f';



%Weight = planner.weightp;%weigth_est;
omega_h = Ramsey.weights.omega_h;     %per history
omegab_h = Ramsey.weights.omegab_h;   %per history
%~ omega_h = Ramsey.weights.omega_param_h; #per history
%~ omegab_h = Ramsey.weights.omegab_param_h; #per history
%%

    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %%%%%%%%%%%     variables and params     %%%%%%%%%%%%
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

str = ['var\n'] ; fprintf(fid, str);

str = ['u r w FK FL A B K mu TFP Ltot Ctot Y\n'] ; fprintf(fid, str);
str = ['ut Kt Ct Bt wt rt taut tauk mut Zt Yt Lt tkt kappa BoYt G l\n'] ; fprintf(fid, str);

%str = ['r w   mu   \n'] ; fprintf(fid, str);
%str = ['tau     \n'] ; fprintf(fid, str);


for h = 1:Ntot
    str = ['a',num2str(h),' '] ; fprintf(fid,str);     %  end-of-period
    str = ['at',num2str(h),' '] ; fprintf(fid, str);   % beginning-of-period
   % str = ['l',num2str(h),' '] ; fprintf(fid, str);   % beginning-of-period
    str = ['x',num2str(h),' '] ; fprintf(fid, str);
    str = ['psib',num2str(h),' '] ; fprintf(fid, str);
    str = ['lambdacb',num2str(h),' '] ; fprintf(fid, str);
    %~ str = ['lambdac'e,num2str(h),' '] ; fprintf(fid, str);
    %~ str = ['lambdal',num2str(p),' '] ; fprintf(fid, str);
    str = ['lambdacbt',num2str(h),' '] ; fprintf(fid, str);    % last period
    %str = ['lambdalb',num2str(h),' '] ; fprintf(fid, str);    % last period
    if mod(h,10)==0; %multiple de 10
        str = ['\n'] ;         fprintf(fid, str);
    end;

end;

str = ['; \n\n'] ; fprintf(fid, str);
str = ['varexo eps; \n\n'] ; fprintf(fid, str);
str = ['parameters\n'] ; fprintf(fid, str);
str = ['beta alpha phi abar delta gamma chi rho_u Gss;\n\n'] ; fprintf(fid, str);

    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %%%%%%%%%%%        initialisation        %%%%%%%%%%%%
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


str = ['alpha','   = ',num2str(alpha,formatSpec),';\n']; fprintf(fid, str);
str = ['beta','    = ',num2str(eco.beta,formatSpec),';\n']; fprintf(fid, str);
str = ['phi','     = ',num2str(eco.phi,formatSpec),';\n']; fprintf(fid, str);
str = ['abar','    = ',num2str(eco.a_min,formatSpec),';\n']; fprintf(fid, str);
str = ['delta','   = ',num2str(eco.delta,formatSpec),';\n']; fprintf(fid, str);
str = ['gamma','   = ',num2str(eco.sigma,formatSpec),';\n']; fprintf(fid, str);
str = ['rho_u','   = ',num2str(rho_u,formatSpec),';\n']; fprintf(fid, str); %HERE
str = ['Gss','     = ',num2str(solution.G,formatSpec),';\n']; fprintf(fid, str);
str = ['chi','     = ',num2str(eco.chi,formatSpec),';\n']; fprintf(fid, str);

%%
equa = 1; % Numbering equations
tol = 0*10^-8;  % threshold for transition, to avoid to consider very small transitions

str = ['\n\n'] ; fprintf(fid, str);

    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %%%%%%%%%%%       model definition       %%%%%%%%%%%%
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


str = ['model;\n\n']; fprintf(fid, str);
for h = 1:Ntot
    
    %%%%%%%% Budget constraint
    str = ['x',num2str(h),' = (1/(chi*taut))*l^(1+1/phi)*(',num2str(trunc.y0_h(h),formatSpec),')^taut + (1+r)*','at',num2str(h),'-a', num2str(h), ...
           '; %% equation ', num2str(equa),'\n']; fprintf(fid, str); equa = equa+1;

    %%%%%%%% Definition of atilde (at)
    str = ['at',num2str(h),' = 0 ']; fprintf(fid, str);
    for hi = 1:Ntot
        if (Pi_h(h,hi))>tol   % hi to h
          str = ['+ ',num2str(S_h(hi)*Pi_h(h,hi)/S_h(h),formatSpec),'*a',num2str(hi),'(-1)']; fprintf(fid, str);
        end
    end
    str = ['; %% equation ', num2str(equa),'\n']; fprintf(fid, str); equa = equa+1;

    %%%%%%%% Euler equation
    if P(h)==1 % not constrained
        str = [num2str(xsis.xsiuE(h),formatSpec),'*(x',num2str(h),')^-gamma =', num2str(trunc.resid_E_h(h),formatSpec),' + beta*(1+r(+1))*(']; fprintf(fid, str);
        for hp = 1:Ntot
            if (Pi_h(hp,h)*xsis.xsiuE(hp))>tol
                str = ['+ ',num2str(Pi_h(hp,h)*xsis.xsiuE(hp),formatSpec),'*(x',num2str(hp),'(+1) )^-gamma']; fprintf(fid, str);
            end
        end
        str = ['); %% equation ', num2str(equa),'\n']; fprintf(fid, str); equa = equa+1;
    else
        str = ['a',num2str(h),' =  ',num2str(trunc.a_end_h(h),formatSpec),';']; fprintf(fid, str);str = [' %% equation ', num2str(equa),'\n']; fprintf(fid, str); equa = equa+1;
    end

    
    %%%%%%%% definition psi_hat_bar
     str = ['psib',num2str(h),' = - mu*',num2str(S_h(h),formatSpec ) ,' + ',num2str(omegab_h(h)*xsis.xsiuE(h),formatSpec ),'*(x',num2str(h),')^-gamma - (lambdacb',num2str(h),...
            ' - (1+r)*lambdacbt',num2str(h), ')*',num2str(xsis.xsiu2(h),formatSpec ),'*(-gamma*(x',num2str(h),')^(-gamma-1)); %% equation ', num2str(equa),'\n']; fprintf(fid, str); equa = equa+1;

    %%%%%%%% FOC individual savings (at): Euler on psib
    %%%%%%%% remark: psib(h) = S(h)*psi(h) = S(h)*sum_hp(psib(hp)/S(hp))
    if P(h)==1
        str = ['psib',num2str(h),'  = ', num2str(S_h(h),formatSpec),'*beta*(1+r(+1))*(0']; fprintf(fid, str);
        for hp = 1:Ntot
            if Pi_h(hp,h)>tol   %be careful transpose from h to hp
                str = ['+ ',num2str(Pi_h(hp,h),formatSpec),'*psib',num2str(hp),'(+1)/',num2str(S_h(hp),formatSpec)]; fprintf(fid, str);
            end
        end
        str = ['); %% equation ', num2str(equa),'\n']; fprintf(fid, str); equa = equa+1;
    else
        str = ['lambdacb',num2str(h),' =  0; %% equation ', num2str(equa),'\n']; fprintf(fid, str); equa = equa+1;
    end

    %%%%%%%% lambdactc : lambdac previous period
    str = ['lambdacbt',num2str(h),' = 0 ']; fprintf(fid, str);
    for hi = 1:Ntot
        if Pi_h(h,hi)>tol  %be careful transpose
            str = ['+ ',num2str(Pi_h(h,hi),formatSpec),'*lambdacb',num2str(hi),'(-1)']; fprintf(fid, str);
        end
    end
    str = ['; %% equation ', num2str(equa),'\n']; fprintf(fid, str); equa = equa+1;
  
end

str = ['\n\n %% FOC planner \n\n']; fprintf(fid, str);



%%%%%%%% FOC planner l
str = ['0 = ((1 +1/phi)/(chi*taut))*(l^(1/phi))*( ']; fprintf(fid, str);
for h = 1:Ntot
    if S_h(h)>tol   %be careful transpose from h to hp
        str = ['+ psib',num2str(h),'*(',num2str(trunc.y0_h(h),formatSpec),')^taut']; fprintf(fid, str);
        if mod(h,10)==0;str = ['\n'] ;fprintf(fid, str); end
    end
end
str = [') - mu*( ']; fprintf(fid, str);
for h = 1:Ntot
    if S_h(h)>tol   %be careful transpose from h to hp
        str = ['+',num2str(S_h(h),formatSpec),'*(l^(1/phi)/chi)*',num2str(trunc.y0_h(h),formatSpec),'^taut',...
               ' -',num2str(S_h(h),formatSpec),'*FL*',num2str(trunc.y0_h(h),formatSpec),'^(1+taut/(1/phi+1))']; fprintf(fid, str);
        if mod(h,10)==0;str = ['\n'] ;fprintf(fid, str); end
    end
end

str = ['); %% equation ', num2str(equa),'\n']; fprintf(fid, str); equa = equa+1;


%%%%%%%% FOC planner B
str = ['mu = beta*(mu(+1)*(1+FK(+1)));  %% equation ', num2str(equa),'\n']; fprintf(fid, str); equa = equa+1;


%%%%%%%% FOC planner r
str = ['0 = 0 ']; fprintf(fid, str);
for h = 1:Ntot
    if S_h(h)>0
        str = ['+ at',num2str(h),'*psib',num2str(h),' + lambdacbt',num2str(h),'*',num2str(xsis.xsiuE(h),formatSpec),'*(x',num2str(h),'^-gamma)']; fprintf(fid, str);
    end
    if mod(h,10)==0;str = ['\n'] ;fprintf(fid, str); end
end
str = [';  %% equation ', num2str(equa),'\n']; fprintf(fid, str); equa = equa+1;


%%%%%%%% FOC planner taut
str = ['\n %% FOC planner tau :', num2str(equa),'\n']; fprintf(fid, str);
str = ['0 = ( (l^(1+1/phi))/(chi*taut) )*(']; fprintf(fid, str);
for h = 1:Ntot
    str = ['+ psib',num2str(h),'*(',num2str(trunc.y0_h(h),formatSpec),'^taut)*(-1/taut + log(',num2str(trunc.y0_h(h),formatSpec),'))'];fprintf(fid, str);
    if mod(h,10)==0;str = ['\n'] ;fprintf(fid, str); end
end
   str = [' )- mu*(l/(1/phi+1))*( '] ;fprintf(fid, str);
for h = 1:Ntot
    str = ['+',num2str(S_h(h),formatSpec),'*log(',num2str(trunc.y0_h(h),formatSpec),')*(  (l^(1/phi))/(chi)*',num2str(trunc.y0_h(h),formatSpec),'^taut'];fprintf(fid, str);   
    str = [' - FL*',num2str(trunc.y0_h(h),formatSpec),'^(1+taut/(1/phi+1))             )'];fprintf(fid, str);       
    if mod(h,10)==0;str = ['\n'] ;fprintf(fid, str); end
end
str = ['); %% equation ', num2str(equa),'\n']; fprintf(fid, str); equa = equa+1;


%%%%%%%%%%%%%%%% Aggregate Constraints %%%%%%%%%%%%%%%%

%%%%%%%% Resource contraint
str = ['G = Gss*(1+u);\n'] ; fprintf(fid, str); str = [' %% equation ', num2str(equa),'\n']; fprintf(fid, str); equa = equa+1;
str = ['G + Ctot + K = Y - delta*K(-1) + K(-1);\n'] ; fprintf(fid, str); str = [' %% equation ', num2str(equa),'\n']; fprintf(fid, str); equa = equa+1;

%str = ['G + (1+r)*B(-1) + (Ctot +A - (1+r)*A(-1))  = Y +B - (r+delta)* K(-1);\n'] ; fprintf(fid, str); str = [' %% equation ', num2str(equa),'\n']; fprintf(fid, str); equa = equa+1;

%%%%%%%% Financial market clearing
str = ['A = B + K; %% equation ', num2str(equa),'\n']; fprintf(fid, str); equa = equa+1;


%%%%%%%% Production side
str = ['FL = (1 - alpha)*TFP*(K(-1)/Ltot)^alpha; %% equation ', num2str(equa),'\n']; fprintf(fid, str); equa = equa+1;
str = ['FK = alpha*TFP*(K(-1)/Ltot)^(alpha-1)-delta; %% equation ', num2str(equa),'\n']; fprintf(fid, str); equa = equa+1;
str = ['u = rho_u*u(-1) + eps;  %% equation ', num2str(equa),'\n']; fprintf(fid, str); equa = equa+1;
str = ['TFP = 1; %% equation ', num2str(equa),'\n']; fprintf(fid, str); equa = equa+1;
str = ['Y = TFP*K(-1)^alpha*Ltot^(1-alpha);  %% equation ', num2str(equa),'\n']; fprintf(fid, str); equa = equa+1;


%%%%%%%% Aggregate labor supply
str = ['Ltot = l*('] ; fprintf(fid, str);
for h = 1:Ntot
    str = ['+', num2str(S_h(h),formatSpec),'*',num2str(trunc.y0_h(h),formatSpec),'^(1+ taut/(1/phi+1))'] ; fprintf(fid, str);
    if mod(h,10)==0;str = ['\n'] ;fprintf(fid, str); end
end;
str = [');   %% equation ', num2str(equa),'\n']; fprintf(fid, str); equa = equa+1;


%%%%%%%% Aggregate consumption
str = ['Ctot = 0'] ; fprintf(fid, str);
for h = 1:Ntot
    str = ['+',num2str(S_h(h),formatSpec),'* (x',num2str(h),'+(( l^(1 + 1/phi)) /(chi*(1/phi+1)))*', num2str(trunc.y0_h(h),formatSpec),'^taut )'] ; fprintf(fid, str);
    if mod(h,10)==0;str = ['\n'] ;fprintf(fid, str); end
end;
str = [';  %% equation ', num2str(equa),'\n']; fprintf(fid, str); equa = equa+1;

%%%%%%%% Wage
str = ['w =  (1/chi)*(1/taut +1/(1/phi+1) ) *l^( (1+1/phi)*(1+1/phi)/(1/phi+1 +taut) )  ;  %% equation ', num2str(equa),'\n']; fprintf(fid, str); equa = equa+1;

%%%%%%%% Taxes
str = ['tauk = 1-r/FK;  %% equation ', num2str(equa),'\n']; fprintf(fid, str); equa = equa+1;
str = ['kappa = w/(FL^( (1/phi+1)*taut/(1/phi+1+taut) ));  %% equation ', num2str(equa),'\n']; fprintf(fid, str); equa = equa+1;


%%%%%%%% Aggregate savings
str = ['A = 0'] ; fprintf(fid, str);
for h = 1:Ntot
    str = ['+', num2str(S_h(h),formatSpec),'*a',num2str(h)] ; fprintf(fid, str);
end;
str = [';  %% equation ', num2str(equa),'\n']; fprintf(fid, str); equa = equa+1;


%%%%%%%%%%%%%%%% log-deviation variables %%%%%%%%%%%%%%%%
% There are denoted with 't': for variable x, xt = (x -x_ss)/x_ss, where x_ss is the steady state value of x.

str = ['ut = 100*u;  %% equation ', num2str(equa),'\n']; fprintf(fid, str); equa = equa+1;
str = ['Kt = 100*(K/',num2str(solution.K(1),formatSpec),'-1);  %% equation ', num2str(equa),'\n']; fprintf(fid, str); equa = equa+1;
str = ['Ct = 100*(Ctot/',num2str(solution.C(1),formatSpec),'-1);  %% equation ', num2str(equa),'\n']; fprintf(fid, str); equa = equa+1;
str = ['Bt = 100*(B/',num2str(solution.B(1),formatSpec),'-1);  %% equation ', num2str(equa),'\n']; fprintf(fid, str); equa = equa+1;
str = ['wt = 100*(w/',num2str(solution.w(1),formatSpec),'-1);  %% equation ', num2str(equa),'\n']; fprintf(fid, str); equa = equa+1;
str = ['rt = 100*(r - ',num2str(solution.R(1)-1,formatSpec),');  %% equation ', num2str(equa),'\n']; fprintf(fid, str); equa = equa+1;
str = ['tkt = 100*(tauk -',num2str(eco.tauk,formatSpec),');  %% equation ', num2str(equa),'\n']; fprintf(fid, str); equa = equa+1;
str = ['mut = 100*(mu/',num2str(Ramsey.lagrangeMult.mu,formatSpec),'-1);  %% equation ', num2str(equa),'\n']; fprintf(fid, str); equa = equa+1;
str = ['Zt = -ut;  %% equation ', num2str(equa),'\n']; fprintf(fid, str); equa = equa+1;
str = ['Lt = 100*(Ltot/',num2str(solution.L,formatSpec),'-1);  %% equation ', num2str(equa),'\n']; fprintf(fid, str); equa = equa+1;
str = ['Yt = 100*(Y/',num2str(solution.Y(1),formatSpec),'-1);  %% equation ', num2str(equa),'\n']; fprintf(fid, str); equa = equa+1;
str = ['BoYt = 100*(B/Y*',num2str(solution.Y/((1+eco.tauc)*solution.A-solution.K),formatSpec),'-1);  %% equation ', num2str(equa),'\n']; fprintf(fid, str); equa = equa+1;
str = ['end;\n\n'] ; fprintf(fid, str);

str = ['%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%\n\n'] ; fprintf(fid, str);
%%%%%%%%%%%%%%%% end of the model part


    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %%%%%%%%%%%     Steady-state model       %%%%%%%%%%%%
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%



str = ['steady_state_model;\n\n'] ; fprintf(fid, str);
for h = 1:Ntot
    str = ['a',num2str(h),'        = ',num2str(trunc.a_end_h(h),formatSpec),';\n'] ; fprintf(fid, str);
    str = ['at',num2str(h),'       = ',num2str(trunc.a_beg_h(h),formatSpec),';\n'] ; fprintf(fid, str);
    %str = ['c',num2str(h),'        = ',num2str(trunc.c_h(h),formatSpec),';\n'] ; fprintf(fid, str);
    str = ['x',num2str(h),'        = ',num2str(x(h),formatSpec),';\n'] ; fprintf(fid, str);
    
    %Here
    %x',num2str(h),'+ ( l^taut) /(chi*(1/phi+1)))*', num2str(trunc.y0_h(h),formatSpec),'^taut )
    
 
    %str = ['l',num2str(h),'        = ',num2str(trunc.l_h(h),formatSpec),';\n'] ; fprintf(fid, str);
    str = ['lambdacb',num2str(h),' = ',num2str(lagMul.lambdacb(h),formatSpec),';\n'] ; fprintf(fid, str);
    %~ str = ['lambdac',num2str(h),'  = ',num2str(lagMul.lambdac(h),formatSpec),';\n'] ; fprintf(fid, str);
    %str = ['lambdalb',num2str(h),' = ',num2str(lagMul.lambdalb(h),formatSpec),';\n'] ; fprintf(fid, str);
    str = ['lambdacbt',num2str(h),' = ',num2str(lagMul.lambdact(h),formatSpec),';\n'] ; fprintf(fid, str);
    str = ['psib',num2str(h),'     = ',num2str(lagMul.psib(h),formatSpec),';\n'] ; fprintf(fid, str);   % recall it is psi multiplied by the size of the bin
end;


str = ['l     = ', num2str(Ramsey.l,formatSpec),';\n'] ; fprintf(fid, str);
str = ['r     = ', num2str(solution.R-1,formatSpec),';\n'] ; fprintf(fid, str);
str = ['w     = ', num2str(solution.w,formatSpec),';\n'] ; fprintf(fid, str);
str = ['FK    = ', num2str(1/eco.beta-1,formatSpec),';\n'] ; fprintf(fid, str);
str = ['FL    = ', num2str(FL,formatSpec),';\n'] ; fprintf(fid, str);
str = ['K     = ', num2str(solution.K(1),formatSpec),';\n'] ; fprintf(fid, str);
str = ['Ltot  = ', num2str(solution.L(1),formatSpec),';\n'] ; fprintf(fid, str);
str = ['mu    = ', num2str(lagMul.mu,formatSpec),';\n'] ; fprintf(fid, str);
str = ['TFP   = ', num2str(1),';\n'] ; fprintf(fid, str);
str = ['A     = ', num2str(solution.A(1),formatSpec),';\n'] ; fprintf(fid, str);
str = ['B     = ', num2str(solution.A(1)-solution.K(1),formatSpec),';\n'] ; fprintf(fid, str);
str = ['u     = ', num2str(0),';\n'] ; fprintf(fid, str);
str = ['tauk  = ', num2str(eco.tauk,formatSpec),';\n'] ; fprintf(fid, str);
str = ['taut   = ', num2str(taut,formatSpec),';\n'] ; fprintf(fid, str);
str = ['Ctot  = ', num2str(solution.C,formatSpec),';\n'] ; fprintf(fid, str);
str = ['Y     = ', num2str(solution.Y,formatSpec),';\n'] ; fprintf(fid, str);
str = ['kappa = ', num2str(eco.kappa,formatSpec),';\n'] ; fprintf(fid, str);

str = ['ut    = 0;\n'] ; fprintf(fid, str);
str = ['Kt    = 0;\n'] ; fprintf(fid, str);
str = ['Ct    = 0;\n'] ; fprintf(fid, str);
str = ['Bt    = 0;\n'] ; fprintf(fid, str);
str = ['wt    = 0;\n'] ; fprintf(fid, str);
str = ['rt    = 0;\n'] ; fprintf(fid, str);
str = ['Tt    = 0;\n'] ; fprintf(fid, str);
str = ['tkt   = 0;\n'] ; fprintf(fid, str);
str = ['mut   = 0;\n'] ; fprintf(fid, str);
str = ['Zt    = 0;\n'] ; fprintf(fid, str);
str = ['Lt    = 0;\n'] ; fprintf(fid, str);
str = ['Yt    = 0;\n'] ; fprintf(fid, str);
str = ['BoYt  = 0;\n'] ; fprintf(fid, str);

str = ['G     = Gss;\n'] ; fprintf(fid, str);
%~ str = ['TLinc = ( G  + (1+FK)*B  - tauk*FK*A - B)/Y ;\n'] ; fprintf(fid, str);
%~ str = ['TKinc = (tauk*FK*A)/Y ;\n'] ; fprintf(fid, str);
%~ str = ['Epsi  = ', num2str(sum(lagMul.psib),formatSpec),';\n'] ; fprintf(fid, str);
%~ str = ['corre = ', num2str(sum(lagMul.psib .* trunc.y0_h)-sum(lagMul.psib),formatSpec),';\n'] ; fprintf(fid, str);

str = ['end;\n\n'] ; fprintf(fid, str);
%%%%%%%%%%% end of the steady state part

    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %%%%%%%%%%%   Steady-state computation    %%%%%%%%%%%
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

str = ['resid;\n\n'] ; fprintf(fid, str);
str = ['steady(nocheck);\n\n'] ; fprintf(fid, str);
 
str = ['check;\n\n'] ; fprintf(fid, str);



    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %%%%%%%%%%%       Aggregate shock        %%%%%%%%%%%%
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

 str = ['shocks;\n var eps =',num2str(G0,formatSpec),';\nend;\n'] ; fprintf(fid, str);



    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %%%%%%%%%%%         Stoch simul          %%%%%%%%%%%%
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    
%~ str = ['stoch_simul (order=1,irf=200) ut Bt tau kappa tauk Yt Ct Kt Lt TLinc TKinc G Y BoYt mut Epsi corre;'] ; fprintf(fid, str);
%~ str = ['stoch_simul (order=1,irf=200) ut Bt tau kappa tauk Yt Ct Kt Lt G Y BoYt mut;'] ; fprintf(fid, str);


str = ['stoch_simul (order=1,irf=500) ut Bt taut kappa tauk Yt Ct Kt Lt;'] ; fprintf(fid, str);

isOctave && fflush(fid);

    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %%%%%%%%%%%        Launch Dynare         %%%%%%%%%%%%
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

    command = ['dynare ',file, ' noclearall'];
    eval(command);

%name = ['fig',num2str(rho_u),'.mat'];
%save (name,'ut_eps','Bt_eps','tau_eps','kappa_eps','tk_eps','Yt_eps','Ct_eps','Kt_eps','Lt_eps','TLinc_eps','TKinc_eps','G_eps','Y_eps','BoYt_eps','mut_eps','oo_')






